Numerical study of metastable states in Ising spin glasses 
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We study numerically the structure of metastable states in the Sherrington-Kirkpatrick spin 
glass. We find that all non-paramagnetic stationary points of the free energy are organized into 
pairs, consisting in a minimum and a saddle of order one, which coalesce in the thermodynamic 
limit. Within the annealed approximation, the entropic contribution of these states, that is the 
complexity, is compatible with the supersymmetry-breaking calculation performed in [A.J. Bray 
and M.A. Moore, J. Phys. C 13 L469 (1980)]. This result indicates that the supersymmetry is 
spontaneously broken in the Sherrington-Kirkpatrick model. 
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In this last year there has been an outburst of new 
interest in the complexity of disordered systems, and in 
particular of spin-glasses 011111111111113. The 
complexity is the entropic contribution due to the expo- 
nentially large number of metastable states, and in mean- 
field disordered models it can be computed by calculating 
the number of minima of the Thouless- Anderson-Palmer 
(TAP) free energy [9J. Metastable states are crucial for 
the dynamical behaviour of disordered systems, and so 
is the complexity. Most notably, in glasses and some 
spin-glasses 0, |n| the complexity triggers a dynamical 
transition which is not associated to any static transition. 
Issues related to the analytic and numerical calculation 
of the complexity are therefore very relevant. 

At the origin of the recent new studies there has been 
the observation done in 2], that the original calcula- 
tion of the complexity for the Sherrington-Kirkpatrick 
(SK) model [12j, performed by Bray and Moore in 
1980 [H |. breaks an intrinsic symmetry of the problem, 
namely a generalized form of the Becchi-Rouet-Stora- 
Tyutin (BRST) supersymmetry Moreover, the 
supersymmetric complexity computed in d was found 
to be much smaller that the supersymmetry-breaking 
complexity of 0]. This discovery reopened, after more 
than twenty years, the problem of how to compute the 
complexity in the SK model, and, more in general, it 
raised the question of whether the BRST supersymme- 
try is spontaneously broken or not in disordered systems. 

At the theoretical level, the situation is presently 
rather open. After the annealed calculation of Q, the 
supersymmetric (SS) complexity has been computed also 
at the quenched level H § , showing that it is consistent 
with the static solution of the SK model and that it coin- 
cides with the one obtained from the Legendre transform 
method at any level of replica symmetry breaking. 
However, it has been shown in that the SS complex- 
ity is stable only at the ground state free energy, where 
by definition it is zero. Thus, were the supersymmetry 
unbroken, this last result would lead to the conclusion 
that in the SK model the complexity of metastable states 
is in fact zero. On the other hand, the supersymmetry- 



breaking (SSB) complexity (which has only been com- 
puted at the annealed level |13| ) is non-zero and stable 
on a finite range of free energy densities. Moreover, it has 
recently been proved in Q that the SSB complexity de- 
scribes pairs of solutions of the TAP equations which are 
either minima or saddles of order one of the free energy. 
According to Q , the two solutions of each pair coalesce 
in the thermodynamic limit. In conclusion, the SS and 
SSB predictions differ considerably. 

Given the ambiguous theoretical situation, a way to 
discriminate between the two pictures above, and thus 
to understand whether or not the BRST supersymmetry 
is spontaneously broken in the SK model, is to perform 
a numerical experiment. Of course, one may question 
what would be the dynamical role of the marginally un- 
stable states described by the SSB complexity, but we 
do not deal with this problem in the present work. Our 
goal here is just to understand what is the structure of 
TAP solutions in the SK model, and which theoretical 
framework correctly describes this structure. 

Our numerical results show that: i) all solutions of the 
TAP equations are either minima or saddles of order one 
of the free energy; ii) the solutions are organized into 
minimum-saddle pairs, connected along a mode which is 
softer the larger the system size; iii) the free-energy dif- 
ference of the paired solutions decreases with increasing 
system size; iv) the annealed complexity of TAP solutions 
as a function of their free energy density is consistent with 
the SSB complexity of 0] . Our results therefore support 
the idea that the BRST supersymmetry is in fact broken 
in the SK model, and that the SSB complexity of is 
correct within the annealed approximation. 

The TAP equations for the SK model are, 
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where m, are the local magnetizations, with i = 1, . . . , N, 
and the random couplings Jij are drawn from a Gaussian 
distribution of mean zero and variance 1/N. These equa- 
tions are obtained by finding the stationary points of the 
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FIG. 1: Number of negative eigenvalues K vs free energy 
density / of 14 non-paramagnetic solutions of a sample at 
N = 80. The solutions are clearly paired into minimum- 
saddle couples, with very close value of /. 



TAP free energy, 

F TAP {m) = ~Yl Ji i mim j ~~J log2 ~ "T^ 1 
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where g = (1/^V) J^i m i ^ s the self-overlap of a solution. 
We solve numerically equations by using the C rou- 
tine BROYDN of Numerical Recipes 01 starting with 
random initial conditions. We analyze 20429 samples 
at N = 20, 5000 samples at N = 30, 1000 samples at 
N = 40, 347 samples at N — 60 and 77 samples at 
N = 80. All data are for T = 0.2 T c . In order to make 
our search of solutions as exhaustive as possible, we use 
two methods. First, all solutions we find in a sample are 
found at least 5 times. Second, we monitor the Morse 
sum, Yla(~ = •"■) wn ere k a is the number of nega- 
tive eigenvalues of solution a. For N — 20, 30, 40, 60 we 
discard samples where this identity is violated for more 
than ±2 (the fraction of discarded samples is always lower 
than 1%). For N = 80, however, due to the low number 
of samples, we do not enforce this last check. 

The first interesting result is that all solutions we find, 
at any value of N, are either minima, or saddles with 
one negative eigenvalue of the Hessian (the paramag- 
net is a minimum). As an example, the solutions of a 
sample at N = 80 are reported in Fig.l. Moreover, all 
non-paramagnetic solutions are organized into minimum- 
saddle pairs: by making a gradient descent starting from 
a saddle and moving along the negative mode, we reach 
in one direction a very close minimum |l8| . and in the 
opposite direction the paramagnetic state. Thus, TAP 
pairs are not directly connected to each other, but they 
are all connected to the central paramagnetic minimum 
in a star-like structure. The self-overlap q of a sad- 
dle is always smaller than that of the connected non- 
paramagnetic minimum. The important point is that 
the two solutions of each pair have very close free energy 
densities /: by plotting the average free energy density 
difference A/ as a function of N, we see that A/ de- 
creases with N (Fig. 2). A power law fit of the data gives 
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FIG. 2: Upper panel: the average free energy density barrier. 
Lower panel: the average modulus of the smallest eigenvalue. 
Lines are power law fits (see text). 



A/ ~ 1/N 1 - 26 . The paired solutions are connected along 
the saddle negative mode, which becomes the smallest 
positive mode of the nearby minimum. The fact that 
the two solutions have very close free energy suggests 
that this direction is extremely flat. In Fig. 2 we also 
plot the average modulus of the smallest eigenvalue A m i n 
as a function of N. A power law fit of the data gives 
I Amin | ~ 1/-/V 140 . Thus the direction connecting the two 
solutions of a pair is flatter the larger N. The exponent 
1/2 proposed in seem to fit the data less satisfyingly. 

These data suggests that for N — > oo the two solu- 
tions of each pair coalesce, forming a single, marginally 
unstable TAP solution. This conclusion is in qualita- 
tive agreement with the theoretical results of results 
found in connection with the SSB annealed complexity. 
It is then natural to ask whether also the annealed nu- 
merical complexity is compatible with the SSB result of 
[l3j |. The annealed complexity is defined as, 



E(/) = ilogjV(/) 



(2) 



where Af(f) is the number of solutions with free energy 
density /, averaged over the disorder J. The SSB an- 
nealed complexity is a smooth function of /, with a 
maximum at a free energy density / max . Let us define 
Smax = S(/ max ) and S;; ax = |S"(/ max )|. For large N 
the total number of TAP solutions M is given by, 
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where we have expanded £(/) to second order close to 
/max- The probability p(f) of finding a solution with free 
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FIG. 3: The variance a of the TAP solutions, an a function of 
N. Full line: the analytic prediction of the SSB calculation, 
a = l/V-N^w- Dashed line: best fit to N~ 1/2 . Dash-dotted 



line: best fit to N~ 1/3 . 



energy density / is therefore given by, 
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which becomes a (^-distribution for N — > oo. For finite N 
the probability distribution p{f) is a numerically directly 
accessible quantity. In particular, its variance a can be 
easily checked without the need of any binning of the 
data. In Fig. 3 we plot the variance a as a function of N. 
The full line is the prediction given by the SSB complex- 
ity, i.e. a = l/VA^max: with £max B) = 8.9. Numerical 
data are compatible with the SSB prediction. Moreover, 
if we perform a fit to a ~ TV -1 / 2 we find S max = 10.3, 
and if we leave the exponent free, fitting to a ~ iV -7 we 
get 7 = 0.58 (not shown). We also show the best fit to 
a r~j TV -1 / 3 , which is the exponent reported in This 
value, however, does not seem compatible with our data. 

By making a (arbitrary) binning of the data in free 
energy, we can compute the full probability p(f) and 
thus our numerical estimate of the complexity, £(/) — 
log[p(f)Af]/N, where N is the total average number of 
solutions per sample we find numerically. Results are 
shown in Fig. 4 and compared to the theoretical curve 
of the SSB complexity £(/). Even though the numer- 
ics is not excellent (especially for N = 80), the data 
seem compatible with the SSB prediction and clearly in- 
dicate that the annealed complexity is nonzero. How- 
ever, a more quantitative comparison, not dependent on 
the binning, is clearly desirable. To do this we study 
the position of the maximum / max and the value of the 
complexity at this point E max , at various N. In order 
to avoid a measurement depending on the binning, we 
study the integral of p(f), i.e. <?(/) = J Q dxp{x), and 
interpolate <?(/) with a cubic function around its inflec- 
tion point. In this way we have an estimate of both / max 
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FIG. 4: Annealed complexity vs free energy density for var- 
ious N; lines connecting the points are only a guide for the 
eye. Full line: analytic SSB prediction. 
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FIG. 5: Upper panel: / max as a function of 1/AT; the dash- 
dotted line is the value of the equilibrium free energy /o. 
Lower panel: E max , as a function of 1/N; the dashed line 
is a 1/N fit to the data. In both panels, the full line is the 
N — oo value of the SSB prediction. 



and S max , reported in Fig. 5 as a function of 1/N. The 
value of / m ax is basically constant, and it agrees quite 
well with the theoretical SSB value = —0.654. 

On the same scale we report the value of the ground 
state (equilibrium) free energy density fa: our data 
do not seem compatible with the law / max — * fo for 
— > oo, proposed in [8|. The value of S max depends 
more strongly on N, as expected from equation @, i.e. 
££ ax = E~„ + log(7V£ max )/(2A0 + 0(1/N). The ex- 
trapolation for A^ — » oo of our data for E max agrees very 
well with the SSB value E^f x B) = 0.052. 

A last important consistency check between our nu- 
merical data and the SSB calculation concerns the be- 
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FIG. 6: The stability parameter x p vs the self-overlap q, para- 
metrically in the free energy density /. Full line: the analytic 
SSB prediction. Data in the box have f < fo- 



haviour of the overlap q and the stability parameter x p , 
defined as, x p = 1 — 1 — mf) 2 /N. In the thermo- 
dynamic limit the condition x p > must be satisfied in 
order to have a physically acceptable TAP state [19|. In 
fact, it is precisely this last condition which is violated by 
the SS complexity in all points but the equilibrium free 
energy / . Both q and x p depend on the free energy den- 
sity / of the solutions. In Fig. 6 we plot x p as a function 
of q, parametrically in /, at various values of 7Y. We see 
that x p is in fact positive for all solutions in the region 
where £(/) > 0. The data agree quite well with the SSB 
prediction, even in the phase / < /o, where £(/) < 0. 

Our algorithm finds more easily solutions with lower 
free energy, as can be seen by checking the frequency each 
solution is found, and this effect is stronger the larger 
N. Therefore, at N = 80 (where the Morse check is 
not enforced) we are probably missing a number of solu- 
tions with / > / max . We believe this is the reason why 
the values of a and / max at N = 80 in Figs. 3 and 5 
are slightly smaller than expected, irrespective of error 
bars. The same problem arises in the reconstruction of 
E(f;N = 80) for / > / max in Fig.4. The problem of 
finding solutions at finite T of the TAP equations is well 
known, and for this reason numerical investigations in 
the past have been done by using alternative methods: 
T = studies of one spin-flip stable states l20l. naive 
TAP equations or modified TAP equations)^. To 
the best of our knowledge, the present work is the first 
extensive study of the full TAP equations at finite T. 

In summary, our results support the conclusion that 
at the annealed level the complexity of the SK model is 
given by the SSB solution of 13], and therefore that the 
BRST supersymmetry is spontaneously broken in this 
system. Considering the fact that the supersymmetry is 
not broken in the p-spin spherical model, it may be ar- 
gued that systems characterized by full replica symmetry 
breaking (as the SK model), and systems solved by one 



step of replica symmetry breaking (as the p-spin), be- 
long to two different classes also for what concerns the 
spontaneous breaking of the supersymmetry. The ori- 
gin of this different behaviour must lie in the geometric 
structure of the states. In particular, the marginality of 
TAP saddle-minima pairs in the SK model is probably 
the cause of the spontaneous supersymmetry breaking. 
This fact is also likely to be connected to the different 
dynamical behaviour of these two classes of systems: the 
dynamics of the p-spin model is heavily influenced by the 
TAP states, while in the SK model, despite the non-zero 
SSB complexity, the dynamics asymptotically reproduces 
the static results, with no role played by the metastable 
states. The supersymmetry could then be an elegant way 
of distinguishing different dynamical classes. 
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